Quantum System Identification by Bayesian Analysis of Noisy Data: Beyond 

Hamiltonian Tomography 



S. G. Schirmer^Q and D. K. L. Oi 2 {D 

1 Department of Applied Maths and Theoretical Physics, 
University of Cambridge, Cambridge CB3 OWA, United Kingdom 
2 SUPA, Department of Physics, University of Strathclyde, Glasgow G4 ONG, United Kingdom 

(Dated: November 6, 2009) 

We consider how to characterize the dynamics of a quantum system from a restricted set of initial 
states and measurements using Bayesian analysis. Previous work has shown that Hamiltonian 
systems can be well estimated from analysis of noisy data. Here we show how to generalize this 
approach to systems with moderate dephasing in the eigenbasis of the Hamiltonian. We illustrate 
the process for a range of three-level quantum systems. The results suggest that the Bayesian 
estimation of the frequencies and dephasing rates is generally highly accurate and the main source 
of errors are errors in the reconstructed Hamiltonian basis. 
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I. INTRODUCTION 

Recent advances in nanofabrication technology in- 
creasingly enable the construction of devices operating 
in the quantum regime. However, to utilize coherence 
effects for practical applications such as quantum infor- 
mation processing and communication tasks requires the 
ability to engineer their dynamics with high precision. 
Considerable progress in the area of laser technology and 
optimal control has shown that precise coherent manip- 
ulation of the dynamics is not infeasible for a variety 
of quantum systems, but such control requires accurate 
knowledge of the system's dynamical behavior and re- 
sponse to external fields, which can be used to construct 
accurate models from which effective control designs can 
be engineered. The problem is particularly acute for 
manufactured systems, due to inevitable variations in the 
manufacturing processes, which ensure that the exact be- 
havior of each device is unique and must be individually 
measured and characterized. 

For the manufacture of large-scale practical devices, 
the design and operation of each device should be as 
simple as possible meaning that the physical resources 
available to initialize and measure the state of a sys- 
tem are usually restricted to a single basis set defined 
by static electrode geometry. In normal operation, any 
state can be produced from an initial fiducial one by ap- 
plying a suitable unitary rotation. This also enables us to 
effectively perform measurements in an arbitrary basis, 
and given both these abilities, one can perform quan- 
tum process tomography [TJ [SJ. However, the problem 
of characterizing a device is not trivial since initially, if 
one does not yet know the control response of the system, 
one cannot generate the unitary rotations required in the 
first place, leading to a Catch-22 situation. 
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What is required is a method of bootstrapping the con- 
trol and characterization process so that the system dy- 
namics and response can be incrementally assessed until 
full control and process tomography is possible, and only 
using the in situ resources. Hence, we have developed 
techniques based upon the analysis of generalized coher- 
ent oscillation data from Rabi or Ramsey-type experi- 
ments. There are several approaches to the analysis of 
such experimental data including frequency-domain and 
time-domain analysis. In the regime of a single system 
transition, Fourier analysis is effective but in the presence 
of multiple signals, it ceases to an optimal estimator. 

In previous work, we have shown how Bayesian signal 
analysis can be effective in determining accurate model 
parameters in generic two-qubit Hamiltonian systems 
where multiple frequencies are present. In this work, 
we extend the technique to systems with dephasing and 
use Bayesian signal analysis to reconstruct the underly- 
ing dynamics, which are now non-unitary. We apply this 
technique to three-level (qutrit) systems and analyze its 
performance for a range of dephasing rates and find that 
as long as coherent dynamics dominate, which would be 
the case for quantum information purposes, signal pa- 
rameters can generally be reliably extracted and the sys- 
tem effectively reconstructed. 

II. OPEN QUANTUM SYSTEMS 

The evolution of a closed quantum system is governed 
by a time-dependent unitary operator U{t) obeying the 
Schrodinger equation. The evolution of an open quan- 
tum system can be highly complicated but under certain 
conditions it can be described by a master equation 

p(t) = -i[H,p(t)] + L D p(t), (1) 

where [A, B] — AB — BA is the usual matrix commu- 
tator, Lfl is a super-operator describing the interaction 
with the environment, and p is a unit-trace positive op- 
erator p on H representing the state of the system. This 
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form of master equation is generally applicable to sys- 
tems interacting with a memory-less (Markovian) reser- 
voir such as an effectively infinite bath, where it can 
be shown that the superoperator Lo(p) takes the form 
Ld(p) — Sfc^I^fe]/ 9 ' where Vjt are operators on H and 
the superoperators T>[Vk] are defined by 



T>[V k ]p = V kP Vl 



pvlv k ). 



(2) 



Under certain conditions we can make further simpli- 
fying assumptions. For example, dissipative effects in 
open systems weakly coupled to an environment are of- 
ten dominated by a certain types of decoherence such as 
pure phase relaxation or population relaxation processes 
such as the spontaneous emission of photons or phonons. 
These types of processes can be described by relatively 
simple master equations. In the case of pure dephasing 
the dissipation superoperator is often determined by a 
single Hermitian operator V. In this case, it is easy to 
show that the master equation simplifies 
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p(t) = -i[H,p(t)]--[V,[V,p(t)}] 



(3) 
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FIG. 1: Simulated measurement traces of a Qutrit. Ideal 
signal trace with no projection noise (top) and with 1000 rep- 
etition samples per time point (bottom). The signal consists 
of three damped sinusoids and represents the probability of 
measuring the system to be in the computational state |0) if it 
was original initialized in |0). At long times, noise dominates 
the signal which leads to an optimal total sampling time. 



Even with these simplifying assumptions on the open sys- 
tem dynamics we see that full system identification now 
requires the identification two generally independent Her- 
mitian operators H and V, which in general means the 
identification of 2(N 2 — 1) real parameters. Fortunately, 
dephasing often acts in the eigenbasis of the Hamiltonian, 
in which case H and V commute and are simultaneously 
diagonalizable, i.e., there exists a basis {le^)} such that 
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where A„ and 7^ are real, and in this case the iden- 
tification problem reduces to finding a joint eigenbasis 
{je^)} and the corresponding eigenvalues A„ and 7„ of 
H and V, respectively. This simplifies the problem. If 
H = diag(A^), V = diag^) and p is the representation 
of the state in a joint eigenbasis of H and V then it is 
easy to see that the master equation ([I]) gives 



dt 



where w^ v = A M - A„ and = 5(7^ - 7^) 2 , i.e. 
have 
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and if W is the unitary basis transformation that maps 
the measurement basis to the joint eigenbasis of H and V , 
then the evolution of the density operator p with respect 
to the measurement basis is given by p(t) = W^p(t)W. 
Thus the evolution is determined by the transition fre- 
quencies oj^v, dephasing rates T ^ and the relation be- 
tween the system and measurement basis W, which are 
to be determined. 



III. EXPERIMENTAL IDENTIFICATION 
PROTOCOL 

As in previous work H G3 EO El 09 03 we as- 
sume that we can prepare and measure the system in 
a fixed set of (orthonormal) computational basis states 
{|1), 1 2), . . . , \N)}, where N is the Hilbert space dimen- 
sion. No other measurements or resources such as non- 
basis states are assumed to be available initially. The 
basic protocol is to prepare the system in a computa- 
tional state, let it evolve for a period of time, then mea- 
sure the probabilities that the system ends up in one of 
the computational basis states, repeating it for different 
times and all computational basis states. The experi- 
mental data thus consists of N 2 time traces, pe,k(t), with 
k,£ = 1,2,3, which represents the probability that the 
system, initially in state \k) is measured in state \£) after 
evolving under the system Hamiltonian for time t. 

When we include dephasing in the Hamiltonian eigen- 
basis, it can be shown that the observable probabilities 
are 



Pk£ 



where the coefficients are 



s ki;u s k£iii C0S(A^ ;A(I/ ), 



Cki = Yju ■ 



'ktv 



(7) 
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(8b) 
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Here s k t,v and 8 k i. v are the amplitude and phase of the 
complex number {l\£ v ) {£ v \k) and A M .^ = 5 M , V - 5 kl .^ 
is the phase difference. 
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FIG. 2: Qutrit Power Spectra. Exponential damping broad- 
ens the peaks leading to difficulties in accurately determining 
their frequencies (green dot-dashed curve). In extremis peaks 
merge (red solid curve) and fewer than three frequencies can 
be observed (blue dashed curve). However the power spec- 
trum can be used as an initial starting point for Bayesian 
estimation. 



If the Hamiltonian is known to be real-symmetric in 
the computational basis, which is the case for many sys- 
tems including atomic and molecular systems, where the 
off-diagonal elements of the Hamiltonian are usually real 
transition strengths or dipolc moments, and spin sys- 
tems, the problem can be simplified. The eigenvectors 
of a real-symmetric matrix are real, thus the phases 5m-,v 
must be multiples of 7r so that e l5fcf; " = ±1, and since 
the sine of a multiple of 7r vanishes, we have b^.^ = 0. 
In many cases the signs of the off-diagonal matrix el- 
ements are also known, e.g. for a spin chain in an 
anti-ferromagnetic material, the off-diagonal elements are 
positive, as the case for many atomic or molecular sys- 
tems. We then have 
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(t) = c M {t) + 2 a M ,^e tr *" cos(uvt) (9) 



with aui-^v = sm-,uSm;ij. and c k g = Y,v s kt,w which fur- 
ther simplifies the reconstruction. 



IV. BAYESIAN PARAMETER ESTIMATION 

This shows that the identification problem for dephas- 
ing that acts in the system's natural basis is similar to 
the Hamiltonian identification problem except that we 
also have to determine the dephasing rates r M „. From 
the measurement results obtained from time traces like 
in Fig. [I] we must extract signal frequencies cu^ and 
damping rates as well as the amplitudes Chi, aki-^ 
and bki-^v in order to be able to perform reconstruction 
of the system dynamics. 



We can do this again by Bayesian estimation, max- 
imizing the likelihood that a particular process gener- 
ated the observed signal [3]. For convenience we la- 
bel the transition frequencies of the system to m , assum- 
ing Lo m +i > w m > 0, and the corresponding dephasing 
rates T m , and define the vectors u> = (u> m ), T = (T m ), 
B-kt = {au-,m) and b w = (b k e;m) where fc,£ range from 1 
to N and m from 1 to the number of transition frequen- 
cies M. According to Eq. [7J the traces should be linear 
combinations 

3 

k£,m 92771 

(t) + c ki (10) 

m=l 

of the mi, = 2M + 1 basis functions 

92m- i(t) = eT tVm cos(w m t), 
52m (t) = e~ tFm s\n(u m t), 



or in the case where H is real-symmetric, the mj, 
basis functions 



(11a) 
(lib) 
(11c) 
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(12a) 
(12b) 



and our objective is to find parameters T m , u> m , akt- m , 
bki-,m and Ck£ that maximize the likelihood of the mea- 
sured data 
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L(a. k t,bu,Ck£,u),a) = J][ n 
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M ex P 



\\Pki 



(13) 

We can eliminate the explicit dependence on the lin- 
ear coefficients a^g, h^g, c k g and the noise variances <Jk£ 
by integration over suitable priors to obtain an explicit 
expression for the probability of a particular model given 
the observed data dkg that depends only on the M tran- 
sition frequencies uj m and corresponding dephasing rates 
T m . Following standard Bayesian analysis [10 we obtain 
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where the averages are defined by 
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(15b) 
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The components hkg- m ar e essentially the orthogonal pro- 
jections of the data onto a set of orthonormal basis vec- 
tors H m (t n ) 
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Ikt:, 



^ H m {t n )dkt-r 



(16) 



71=1 



4 



derived from the (non-orthogonal) basis functions g m {t) 
defined above, evaluated at the respective sample times 



*„., via 



Hm (^n) 
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(17) 



where e m ' m is a rrif, x mi, matrix whose columns e m are 
the normalized eigenvectors — Ge m — a m e m — of the 
?7if, x mi, matrix G — (G mim2 ) with 
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N 

in im 2 ^ ' Hi 
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(18) 



Thus, the parameter estimation problem for a system 
with decoherence acting in the Hamiltonian basis is sim- 
ilar to that for a Hamiltonian system, except that the 
sine and cosine basis functions for the Bayesian analy- 
sis must be modified to damped sinusoids with unknown 
damping rates. 

The objective is to find the frequencies u) and damping 
rates T that maximize -P(u>, r|d^), or equivalently, the 
log-likelihood function 



log 10 P(u;,r|d fc 



mb — N 



ki ) 



N 

£ 

k.e=i 



logi 
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(19) 

Given a solution u> and T that maximizes this log- 
likelihood, it can be shown that the corresponding op- 
timal coefficients in the general case (11) are 



&kt = {{xu;i), (xke-3), ■ ■ ■ , (xk£;m B -2)) , (20a) 
b k t = ({Xkl;2), {xm-a), {Xki;m B -l)) i ( 20b ) 
Ckl = (Xki;m B )l ( 20c ) 

where (£fc^ ;m ) is shorthand notation for the expecta- 
tion values E(xki: m \ijJ, I\ dki) of the linear coefficients 
of the basis functions, given the optimal frequencies a; 
and damping rates T and the data d^e . Similarly in the 
special case (fl2l 
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(21a) 
(21b) 



Since the log-likelihood function is sharply peaked 
with generally many local extrema, finding the global 
optimum using gradient-type optimization algorithms 
starting with a completely random guess for uj and 
r is inefficient. A global optimization such as pat- 
tern search or evolutionary algorithms might circum- 
vent this problem, but neither proved either very effec- 
tive in our case, especially for higher-dimensional search 
spaces. Alternatively, starting with a somewhat reason- 
able initial guess, especially for the frequencies, a stan- 
dard quasi-Newton optimization method with cubic line 
search [TTJ H3J E] proved generally very effective in 
finding the global maximum. 



To obtain an initial estimate for the frequencies we 
used the sum of the power spectra of the signals. Al- 
though the peaks in the power spectrum are not optimal 
frequency estimators when there are multiple frequencies 
and the exact peak locations can be difficult to ascertain 
even for systems with only three frequencies, as Fig. [2] 
shows, rough estimates of the peak locations usually seem 
to provide a reasonable initial guess for the gradient- 
based likelihood optimization routine. In principle the 
damping rates could be estimated from the peaks widths 
as well but these estimates can be tricky, especially for 
overlapping and minor peaks, hence we chose multiple 
runs with random initial guesses for the damping rates 
r and selected the run with the highest final likelihood 
( "global" maximum) . 

Given the extracted signal parameters we have to solve 
two further inverse problems: (i) reconstructing the level 
structure from the frequencies and (ii) constructing the 
matrix W that relates the Hamiltonian basis to the com- 
putational basis. The former usually involves analyzing 
the relationships between the frequencies as illustrated 
in [S]. In general this is be tricky but for a qutrit system, 
is analysis is essentially trivial. The basis reconstruction 
requires solving further optimization problems to find the 
coefficients sut,^ such that Eqs (JsJ> are satisfied given the 
estimates for the parameters ay ; „, b^-a, Cki and A^.^ 
derived in the previous step. Due to finite sampling and 
noise, the inversion may not be exact, hence we recast 
it as a constrained optimization problem and solve it as 
described in [S]. 

Our previous analysis [S] also shows that we can only 
identify a single generic Hamiltonian up to equivalence 



H ~ D^HD + XI, 



(22) 



where D = diag(l, e lSl2 , . . . , e tSlN ) is a diagonal unitary 
matrix, in the basis of the measurement. However, if 
the off-diagonal elements in the Hamiltonian are known 
to be real and positive, for instance, then the Hamilto- 
nian will be uniquely determined up to a global energy 
level shift Al, at least in the generic case, since we have 
\Hf.g\ = \Hki\ f° r k £. For a quantum control situa- 
tion, the system dynamics can be controlled and hence 
different Hamiltonians can be applied, and in the case 
subsequent Hamiltonians can be fully determined up to 
the gauge fixed by the initial Hamiltonian. By varying 
control parameters and tracking the change in the system 
dynamics, a dynamical control model can be built of the 
system. 



V. RESULTS 

We randomly generated 100 real-symmetric qutrit 
Hamiltonians and dephasing operators with different 
spectral properties and the geometric average of the sys- 
tem Q-factors ranging from 12 to 72. From these we 
generated various data traces corresponding to the stro- 



boscopic sampling described in section III We considered 
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TABLE I: Median Likelihoods (L) and Error Rates (e) for 
Qutrit Systems. For the 100 qutrit systems we compared the 
case with and without dephasing (superscript H) for different 
samples (A e ) per data point. With no sampling noise N^, 
there was a small change in the median errors. For the Aiooo 
case, the median errors increase due to the sampling noise, 
the addition of dephasing increases the final error by an order 
of magnitude to the 1% region. A simple adaptive scheme 
N var does similarly. The Hamiltonian is reconstructed using 
several runs of the optimization routine, and the solution with 
the minimum basis error is chosen. 



three cases, the zero noise case (N^ ^infinite samples 
per point) , fixed finite sampling with N e = A^iooo = 1000 
experimental repetitions per time point, and an adap- 
tive sampling strategy N var which varies the number of 
samples per point to reach an estimated target signal to 
noise ratio of (pkt{t)) > 10/^/Ne for all k,£ and t with 
an upper limit of N e < 10, 000 for each data point. We 
then applied our parameter estimation and reconstruc- 
tion algorithms to the resulting data traces. A range of 
dephasing rates was studied to see the effect on the re- 
construction of the Hamiltonian part of the dynamics. 
For the purposes of control, accurate determination of 
the Hamiltonian is much more important than a precise 
determination of the dephasing rate, usually it suffices to 
know that they are below certain limits. 

Table U shows the median errors for various cases. 
Comparing the dephasing/no dephasing cases, the errors 
arc similar in the absence of projection noise (A^). The 
frequency u> estimation is slightly more accurate but es- 
timation of the signal amplitudes a^i-^ is slightly less 
accurate since the basis functions depend on T, hence 
errors in both u> and T contribute to errors in the co- 
efficients ciki^v For reduced signal to noise, dephas- 
ing decreases the maximum likelihood and increases fre- 
quency, basis and reconstructed Hamiltonian errors with 
a marked increase in median of the amplitude errors. 
Adaptive sampling overall increases the accuracy of the 
parameter estimation step and the reconstructed Hamil- 
tonian for both Hamiltonian and dephasing systems but 
the improvement is more pronounced for dephasing sys- 
tems. This may be due to adaptive sampling being more 
beneficial for small signal amplitudes i.e., decaying sig- 
nals. This suggests the use of adaptive sampling to in- 
crease the signal to noise ratio for samples at increasing 
times. Alternatively, the sample data can be weighted to 
give precedence to earlier samples. Further exploration 
of these methods will be the subject of future study. 
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FIG. 3: Hamiltonian vs Basis Reconstruction Error S for var- 
ious samplings shows a strong correlation and suggests that 
the total error in the Hamiltonian is dominated by errors in 
the basis reconstruction step that comes about from the sep- 
arate optimization of each basis function from the amplitude 
estimation, which may not lead to orthogonal data vectors. 



Dephasing leads to a reduction in signal at long times 
which can lead us to fitting noise. For strong dephas- 
ing, this leads to reduced accuracy in the estimation of 
the frequencies, and hence increased errors in the other 
parameters. The spread in the Fourier peaks can also 
lead to problems for closely spaced frequencies. This in 
itself is not a problem per se for the Bayesian parameter 
estimation step jS] , except that it can lead to inaccurate 
initial search parameters coming from peak detection in 
the power spectrum. This can be obviated somewhat by 
trying different initial parameters assuming that either 
of the two remaining peaks were doublets and using the 
most likely result. 

For systems of interest for quantum information pro- 
cessing, the dephasing rates should be sufficiently low so 
that the damping of the Rabi-type oscillations do not im- 
pact the scheme greatly. For very small dephasing rates, 
However, it can be a problem if the algorithm overesti- 
mates the dephasing rates which means that the basis 
functions used are not suitable, and this is reflected in 
errors of the estimated amplitudes. For such systems, it 
is a simple enough matter to test models which are purely 
Hamiltonian to see which gives the larger likelihood. 

One factor which limits the reconstruction is that we 
may obtain a set of N x N matrices P v — (ske-v), which 
ideally should be projectors onto orthogonal eigenspaces, 
but may not always form an orthogonal set of projectors. 
We can quantify this basis error by 



S 



max|Tr(P !/ t P /J ) 



(23) 



Fig. [3] shows that there is a strong correlation between S 
and the (relative) error in the final reconstructed Hamil- 
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tonian. Thus, we can use S to choose the best recon- 
structed Hamiltonian from multiple optimization runs 
and as a rough indication of the likely accuracy of the 
reconstructed Hamiltonian. The data also suggests that 
there is little direct correlation between the likelihood 
and errors in the parameter estimation step and the fi- 
nal Hamiltonian error, suggesting that the final error in 
the Hamiltonian is dominated by errors in the basis re- 
construction step. The reconstruction step obviously de- 
pends on the parameter estimates obtained in the first 
step, and poor estimates for the parameters will generally 
result in large Hamiltonian errors, but in some cases the 
basis reconstruction produces poor results even when the 
individual errors in the estimated parameters are small. 
It should be possible to improve the reconstruction step 
by solving the N optimization problems for the Skt- ttl 
simultaneously rather than independently and enforcing 
orthonormality constraints for the basis vectors, but do- 
ing so would require solving a rather more complicated 
optimization problem with several nontrivial constraints. 

VI. DISCUSSION 

Other researchers have also begun to address the prob- 
lem of system characterization with limited resources. 
For example, Leghtas et al. jTS] also consider estimating 
parameters of three-level quantum systems using weak 
continuous population measurements. However, in their 
case it is assumed that most of the system is already 
known including the transition frequencies and the pre- 
cise structure of the Hamiltonian, and there is no intrinsic 
decoherence. They consider extracting only two real pa- 
rameters of the system, the dipole transition strengths 
between levels 1-2 and 2-3, which simplifies the problem 
enormously. 



Burgarth et al. [THl HZ] also consider Hamiltonian char- 
acterization with restricted resources for Heisenberg spin 
chains where only a small subset of spins are individually 
addressable. The form and structure of the Hamiltonian 
is known a priori to be of a particular class, and only the 
coupling strengths and anisotropy of the system Hamil- 
tonian are to be determined. The sign of the couplings 
is also known beforehand. Characterization is achieved 
in this case by preparing different initial states of the 
first spin, letting the system evolve and then performing 
quantum state tomography on the accessible spins. If we 
consider a system of three spins, the first excitation sub- 
space acts as a qutrit. Our protocol could be applied to 
this problem with some modifications. Our scheme does 
not require state tomography, only the determination of 
position of the up-spin, and there is no requirement to 
know the network topology. It would be interesting to ex- 
plore Bayesian analysis of the response of such systems 
for Hamiltonian characterization, and especially the role 
of topology in idcntifiability, and whether it is possible 
to relax the requirement for addressability of all spins. 

In summary, we have shown that our current two-step 
procedure of Bayesian parameter estimation followed by 
a reconstruction via optimization works in the presence of 
dephasing on three-level systems. However, we find that 
the reconstruction step is a weak point of our current im- 
plementation. It may be possible to eliminate the param- 
eter estimation step and directly apply Bayesian maxi- 
mum likelihood estimation upon the dynamical system 
parameters. This would have the advantage of always 
giving admissible solutions at all steps. Another direc- 
tion which should be explored is adaptive sampling, not 
only varying experimental repetitions per data point, but 
also using non-uniform time-domain sampling for better 
frequency discrimination. 
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